home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / svdfit.pro < prev    next >
Text File  |  1997-07-08  |  8KB  |  250 lines
  1. ; $Id: svdfit.pro,v 1.12 1997/02/04 13:11:51 stevep Exp $
  2.  
  3. FUNCTION SVDFIT,X,Y,M, YFIT = yfit, WEIGHTS = weights, CHISQ = chisq, $
  4.     SINGULAR = sing, VARIANCE = variance, COVAR = covar, SIGMA=sigma, $
  5.     Function_name = function_name, A=A,LEGENDRE=legendre,DOUBLE=double
  6. ;+
  7. ; NAME:
  8. ;    SVDFIT
  9. ;
  10. ; PURPOSE:
  11. ;    Perform a general least squares fit with optional error estimates.
  12. ;
  13. ;    This version uses the Numerical Recipies (2nd Edition) function
  14. ;    SVDFIT.  A user-supplied function or a built-in    polynomial or
  15. ;    legendre polynomial is fit to the data.
  16. ;
  17. ; CATEGORY:
  18. ;    Curve fitting.
  19. ;
  20. ; CALLING SEQUENCE:
  21. ;    Result = SVDFIT(X, Y, [M])
  22. ;
  23. ; INPUTS:
  24. ;    X:    A vector representing the independent variable.
  25. ;
  26. ;    Y:    Dependent variable vector.  This vector must be same length 
  27. ;        as X.
  28. ;
  29. ; OPTIONAL INPUTS:
  30. ;
  31. ;    M:    The number of coefficients in the fitting function.  For 
  32. ;        polynomials, M is equal to the degree of the polynomial + 1.
  33. ;               If not specified and the keyword A is set, THEN
  34. ;        M = N_ELEMENTS(A).
  35. ;
  36. ; INPUT KEYWORDS:
  37. ;
  38. ;        A:  The inital estimates of the desired coefficients. If M
  39. ;               is specified, THEN A must be a vector of M elements. 
  40. ;               If A is specified, THEN the input M can be omitted and
  41. ;               M=N_ELEMENTS(A). If not specified, the initial value
  42. ;               of each coefficient is taken to be 1.0. If both M and A are
  43. ;               specified, them must agree as to the number of paramaters.
  44. ;
  45. ;    DOUBLE:    Set this keyword to force double precision computations. This
  46. ;               is helpful in reducing roundoff errors and improves the chances
  47. ;               of function convergence.
  48. ;
  49. ;    WEIGHTS:    A vector of weights for Y[i].  This vector must be the same
  50. ;        length as X and Y.  If this parameter is ommitted, 1's
  51. ;        (No weighting) are assumed.  The error for each term is
  52. ;        weighted by Weight[i] when computing the fit.  Gaussian or
  53. ;        instrumental uncertianties should be weighted as
  54. ;        Weight = 1/Sigma where Sigma is the measurement 
  55. ;               error or standard deviations of Y. For Poisson or statistical 
  56. ;               weighting use Weight=1/Y, since Sigma=sqrt(Y).
  57. ;
  58. ; FUNCTION_NAME:
  59. ;        A string that contains the name of an optional user-supplied 
  60. ;        basis function with M coefficients. If omitted, polynomials
  61. ;        are used.
  62. ;
  63. ;        The function is called: R=SVDFUNCT(X,M)
  64. ;
  65. ;        where X and M are  scalar values, and the function value is an 
  66. ;        M element vector evaluated at X with the M basis functions.  
  67. ;               M is the degree of the polynomial +1 if the basis functions are
  68. ;               polynomials.  For example, see the function SVDFUNCT or SVDLEG,
  69. ;               in the IDL User Library:
  70. ;
  71. ;        For more examples, see Numerical Recipes in C, second Edition,
  72. ;               page 676-681.
  73. ;
  74. ;        The basis function for polynomials, is R[j] = x)^j.
  75. ;
  76. ;            The function must be able to return R as a FLOAT vector or
  77. ;               a DOUBLE vector depending on the input type of X.
  78. ;        
  79. ;     LEGENDRE: Set this keyword to use the IDL function SVDLEG in the lib
  80. ;               directory to fit the data to an M element legendre polynomial.
  81. ;               This keyword overrides the FUNCTION_NAME keyword.
  82. ;
  83. ; OUTPUTS:
  84. ;    SVDFIT returns a vector of the M coefficients fitted to the
  85. ;    supplied function.
  86. ;
  87. ; OPTIONAL OUTPUT PARAMETERS:
  88. ;
  89. ;    CHISQ:    Sum of squared errors multiplied by weights if weights
  90. ;        are specified.
  91. ;
  92. ;    COVAR:    Covariance matrix of the coefficients.
  93. ;
  94. ;    VARIANCE:    Sigma squared in estimate of each coeff(M).
  95. ;               That is sqrt(VARIANCE) equals the 1 sigma deviations
  96. ;               of the returned coefficients.
  97. ;
  98. ;      SIGMA:   The 1-sigma error estimates of the returned parameters,
  99. ;               SIGMA=SQRT(VARIANCE).
  100. ;
  101. ;    SINGULAR:    The number of singular values returned.  This value should
  102. ;        be 0.  If not, the basis functions do not accurately
  103. ;        characterize the data.
  104. ;
  105. ;    YFIT:    Vector of calculated Y's.
  106. ;
  107. ; COMMON BLOCKS:
  108. ;    None.
  109. ;
  110. ; SIDE EFFECTS:
  111. ;    None.
  112. ;
  113. ; MODIFICATION HISTORY:
  114. ;    Adapted from SVDFIT, from the book Numerical Recipes, Press,
  115. ;    et. al., Page 518.
  116. ;    minor error corrected April, 1992 (J.Murthy)
  117. ;
  118. ;    Completely rewritten to use the actual Numerical Recipes routines
  119. ;    of the 2nd Edition (V.2.06). Added the DOUBLE, SIGMA, A, and
  120. ;    LEGENDRE keywords. Also changed Weight to Weights to match the
  121. ;    other fitting routines.
  122. ;
  123. ;       EXAMPLE:
  124. ;-
  125.     ON_ERROR,2        ;RETURN TO CALLER IF ERROR
  126. ;
  127. ;       Ensure that the proper number/type of parameters are used
  128. ;
  129.         np=N_PARAMS()
  130.         HAVE_IN_A=KEYWORD_SET(A)
  131.         IF HAVE_IN_A THEN IN_A=A
  132.         IF not KEYWORD_SET(LEGENDRE) THEN use_legendre=0 else use_legendre=1
  133.  
  134. ;
  135. ;       Either two or three parameters must be supplied
  136. ;
  137.         IF np NE 2 and np NE 3 THEN $
  138.       MESSAGE,'Incorrect number of arguments for SVDFIT,x,y,[m]'
  139.  
  140. ;
  141. ;       IF 2 paramaters are supplied, THEN A must be supplied
  142. ;
  143.         IF np EQ 2 THEN $
  144.           IF N_ELEMENTS(A) eq 0 THEN $
  145.          MESSAGE,'The keyword A must be set to call SVDFIT as SVDFIT,x,y' $
  146.           ELSE M=N_ELEMENTS(A)
  147. ;
  148. ;
  149.     IF np EQ 3 THEN BEGIN
  150.       IF (SIZE(M))[0] NE 0 THEN $
  151.             MESSAGE,'The input M must be a scalar.' 
  152.       IF not KEYWORD_SET(A) THEN A=REPLICATE(1.0,M) else $
  153.          IF N_ELEMENTS(A) NE M or (SIZE(A))[0] NE 1 THEN $ 
  154.             MESSAGE,'The keyword A must be an M element vector.' 
  155.         ENDIF
  156.  
  157.     IF (SIZE(X))[0] NE 1 THEN $
  158.             MESSAGE,'The input X must be a vector.' 
  159.  
  160.     IF (SIZE(Y))[0] NE 1 THEN $
  161.             MESSAGE,'The input Y must be a vector.' 
  162.      
  163. ;
  164. ;       Evaluate/Set/Check a few essential paramaters
  165. ;
  166.  
  167.     THRESH = 1.0E-9        ;Threshold used in editing singular values
  168.  
  169.     ndata = N_ELEMENTS(x)     ;SIZE of X
  170.  
  171.     IF ndata NE N_ELEMENTS(y) THEN  $ 
  172.       message, 'Error: X and Y must have same # of elements.'
  173.  
  174.     IF N_ELEMENTS(weight) NE 0 THEN BEGIN
  175.         IF N_ELEMENTS(weight) NE ndata THEN $
  176.           MESSAGE, 'Error: Weights must have the number of elements as X and Y.'
  177.       sig = 1/weight    ;Apply weights
  178.     ENDIF ELSE sig = FLTARR(ndata) > 1.0    
  179.  
  180.     IF n_elements(FUNCTION_NAME) EQ 0 THEN FUNCTION_NAME='svdfunct' 
  181.     IF FUNCTION_NAME EQ '' THEN FUNCTION_NAME='svdfunct' 
  182.         IF use_legendre THEN FUNCTION_NAME='svdleg'
  183.  
  184. ;
  185. ;       If x or y or sig or a is double precision, set use_double to true
  186. ;
  187.         If (reverse(size(a)))[1] eq 5 or $
  188.            (reverse(size(x)))[1] eq 5 or $
  189.            (reverse(size(y)))[1] eq 5 or $
  190.            (reverse(size(sig)))[1] eq 5  then use_double=1
  191.  
  192.         IF not KEYWORD_SET(DOUBLE) THEN use_double=0 else use_double=1
  193.  
  194.        if use_double then begin
  195.         chisq=0.0D
  196.         xx=double(x) & yy=double(y) & a=double(a)
  197.         ssig=abs(double(sig)) > 1E-12
  198.        endif else begin
  199.         chisq=0.0
  200.         xx=float(x) & yy=float(y) & a=float(a)
  201.         ssig=abs(float(sig)) > 1E-6
  202.        endelse
  203.  
  204. ;
  205. ;       Call the actual NR routine.
  206. ;   
  207. ;       Warning, direct use of this function is not supported, as
  208. ;       the calling sequence is scheduled to change in a future release
  209. ;       of IDL.
  210. ;
  211.         NR__SVDFIT,FUNCTION_NAME,xx,yy,ssig,ndata,A,M,COVAR,CHISQ,$
  212.             DOUBLE=use_double
  213.       
  214. ;
  215. ;       The variance is the diagonal elements of the covarianec matrix
  216. ;
  217.         variance=COVAR[lindgen(M)*(M+1)];
  218.         sigma=sqrt(abs(VARIANCE))
  219.  
  220.         small=where(variance le max(variance)*thresh,cc)
  221.         IF cc GT 0 THEN variance[small]=0
  222.  
  223.     good = where(variance GT 0, ng) ;Cutoff for sing values
  224.     sing = M - ng        ;# of singular values
  225.     IF sing NE 0 THEN BEGIN
  226.         message, 'Warning: ' + strcompress(sing, /REMOVE) + $
  227.             ' singular values found.',/INFORM
  228. ;
  229.         IF ng EQ 0 THEN return,undefined
  230.     ENDIF
  231. ;
  232. ;     calculate YFIT 
  233. ;
  234.           IF KEYWORD_SET(use_double) THEN BEGIN
  235.              yfit=DBLARR(ndata)
  236.              afunc=DBLARR(M)
  237.           ENDIF ELSE BEGIN
  238.              yfit=FLTARR(ndata)
  239.              afunc=FLTARR(M)
  240.           ENDELSE
  241.          FOR i=0,ndata-1 DO $
  242.         YFIT[i]=total(A * call_function(FUNCTION_NAME, xx[i], M))
  243. ;
  244. ;     Return the fitted parameters
  245. ;
  246.          OUT_A=A
  247.          IF HAVE_IN_A THEN A=IN_A
  248.      return,OUT_A
  249. end
  250.